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UNSTEADY LAMINAR BOUNDARY LAYER CALCULATIONS ON OSCILLATING 


CONFIGURATIONS INCLUDING BACKFLOW 
PART I: FLAT PLATE, OSCILLATING IN ITS OWN PLANE 

W. Geissler* 

Ames Research Center 

SUMMARY 


A finite difference method has been developed to calculate the unsteady boundary • 
layer over an oscillating flat plate. Low- and high-frequency approximations were 
used for comparison with numerical results. Special emphasis was placed on the 
behavior of the flow and on the numerical calculation procedure as soon as reversed 
flow has occurred over part of the oscillation cycle. The numerical method displayed 
neither problems nor singular behavior at the beginning of or within the reversed 
flow region. Calculations, however, came to a limit where the back-flow region 
reached the plate's leading edge in the case of high oscillation amplitudes. It is 
assumed that this limit is caused by the special behavior of the flow at the plate's 
leading edge where the boundary layer equations are not valid. 

INTRODUCTION 


Unsteady viscous flows play a dominant role in several fields of aeronautical 
applications and have severe influences on lift and drag characteristics, particularly 
in cases where flow separation occurs. Dynamic-stall investigations on oscillating 
3-D rotor-blade tips (ref. 1) and active control investigations on a supercritical 
wing section with an oscillating, trailing-edge control (ref. 2) are two of these 
cases in which unsteady separated flows are dominant. These two cases have been 
investigated recently by using intensive wind-tunnel measurements made in the Deutsche 
Forschungs- und Versuchsanstalt fUr Luft- und Raumfahrt (DFVLR) GUttingen. Both cases 
have been compared with inviscid-panel-method calculations (refs. 1, 3) to illustrate 
the influence of viscosity and its effects on the steady and unsteady airloads under 
nonseparated and separated flow conditions. These studies indicated that the influ- 
ence of viscosity in nonseparated flow is created mainly by boundary layer displace- 
ment effects. These effects can easily be taken into account in theoretical calcula- 
tions by adding the displacement thickness to the wing profile and then calculating 
pressures and forces on the modified configuration. 

Under unsteady separated flow conditions, however, this procedure is not appli- 
cable. Comparisons between inviscid theory and experiments in this case have shown 
that special characteristic features of the pressure distributions occur in most of 
the investigated cases. These features can be described by strong phase shifts 
within the first-harmonic unsteady pressure distributions when they are compared to 
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inviscid calculations. These phase shifts may lead to a shift of the real-part pres- 
sure peak (at the hinge line of a wing having oscillating control) into the imaginary 
part, exerting a strong influence on the applicability of wing-control systems for 
active-control purposes. 

To investigate these viscous phenomena, unsteady boundary-layer calculations 
have been performed which include flow regimes in which back-flow occurs along the 
oscillating wall. The present study includes the case of a semi-infinite flat plate, 
oscillating in its own plane. This problem has been treated in several investiga- 
tions from the first studies by Lighthill (ref. 4), who calculated approximate high- 
and low-frequency solutions to the more sophisticated numerical procedures of Phillips 
and Ackerberg (ref. 5), which deal with regions of backflow. McCroskey and Philllppe 
(ref. 6) extended their method to the dynamic stall problem of oscillating profiles 
which include laminar and turbulent unsteady boundary layer calculations. The present 
study of an oscillating flat plate is assumed to be a necessary step in checking the 
validity and applicability of the numerical method which will be applied to the more 
complicated case" of oscillating profiles. Some interesting features have been found 
in the flat plate boundary layer studies which may give an indication of the behavior 
of more complicated unsteady viscous flows. 

( 


UNSTEADY BOUNDARY LAYER EQUATIONS 


Unsteady, viscous-incompressible flow in the vicinity of an oscillating wall is 
governed by the well-known set of unsteady boundary layer equations expressing the 
conservation of mass and momentum (ref. 7): 


5u _3v 

ax an 


= 0 


( 1 ) 


8u ' _8u , - _3u _ 3 jd a 2 u 

ai ax an = “ 3x 9ti 2 

with -(3p/3x) = (au/3T) in the flat plate case. 


( 2 ) 


In equations (1) and (2), the different quantities have been made dimensionless, 
and the normal coordinate and the normal velocity have been multiplied by v^Re. The 
boundary conditions for the system of equations read: 

u = v = 0 for n = 0 

u = U for n = 00 


In the flat plate case, the velocity at the outer edge of the boundary layer is 
assumed to be given by 


U(T) = U w (l + e • e il0 * T ) (3) 

with e as the oscillation amplitude. Three different solution procedures will be 
developed to include: 
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1. A low-frequency approximation, where the unsteady parameter E, = x • w* < 1; 

2. A high-frequency approximation, where the unsteady parameter £ » 1; 

3. A finite-difference procedure which numerically calculates the unsteady 
boundary layer for arbitrary values of £. 

Comparisons between the low- /high-frequency approximations and the numerical results 
will provide an excellent check of the validity of the numerical procedure. 

LOW-FREQUENCY APPROXIMATION 


In the boundary layer equations (1) and (2), only the steady and first harmonic 
terms are retained. The stream-wise velocity component is defined by 

u = f’(n) + £$’ (£,n)e iaJ * T (4) 

As long as E, < 1, a power-series expansion in E, can be applied for 4: 

*(€,n) = E eV(n)e iw * T (5) 

k=o k 


leading to 

i 


u = f ' (n) + e 



5 k <J^(u) 


itj*T 

e 


( 6 ) 


with f as the well-known Blasius function. With a properly chosen normal velocity 
component 


--hi !{ 


(nf' - f) + e(n E C <J> 


(n E 

\ k=o 


-k, » 


- E 

k=o 




) 


iu*T 


- 2eg 


E 

k=i 


k? k_1 <f>i 


iu*Tl 


(7) 


the continuity equation (1) is automatically fulfilled. 


The following system of equations is obtained by combining expressions (6) and 
(7) and the corresponding derivatives with equation (3) in the momentum equation (2): 
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Equation 

f • f" + 2f”’ = 0 
(Blasius equation) 

\ «♦; +\ f"*. + ♦!" 

(quasi-steady) 

' i - *5 - + \ 

+ | 

(linear) 

- ♦;_! - nf l *i +i 

+ (n + |) £’’* n + +A" 

(quadratic and higher) 


Boundary condition 


f - f' = 0 

V = 1 

to - K - 0 
*5 - i 

<h = <K ■ o 

<(>I = o 

♦n - < - 0 

♦A - 0 


for n = 0 

for n = 00 

for n = 0 

for n = 08 

for n = 0 

for r) = » 

for n = 0 

for ri = 00 


( 8 ) 


The system of equations (8) is solved sequentially, starting with the Blasius equa- 
tion and then solving the quasi-steady, linear, and higher-order equations by taking 

the solutions of the previous steps into account. 

. ' / 

The numerical solution of the system of ordinary differential equations (8) is 
carried out by a fourth-order Runge-Kutta integration procedure which has been used 
for calculating steady, 3-D stagnation-point flows (refs. 8, 9). Typical results are 
plotted in figure 1, which gives the real and imaginary parts of the function $ 

(eq. (5)) versus n for various parameters 5. For each of these calculations, the 
series expansion in equation (5) was represented up to the 11th member. The quasi- 
steady curve in figure 1 (broken line) is the limiting case of k = 0 in equation (5), 
which indicates that the function $ is independent of 5 and therefore is indepen- 
dent of oi*. The overshoot at the outer boundary condition, which increases with 
increasing 5> is typical behavior of the real parts of , Further results of the 
low-frequency approximation are compared directly with the numerical results. 


HIGH-FREQUENCY APPROXIMATION 


In the high-frequency case, a procedure similar to that used for low frequencies 
is not straightforward. A representation of $ by a series expansion of 1/5, for 
instance, does not lead to the correct approximation of the governing equations. 
Therefore, a different procedure is followed here. 

By using equation (4) for u and using equation (7) for v without taking into 
account the series expansion of 0, the following well-known differential equation 
for $ (ref. 7, p. 434) is derived 

$>”' + ft" - gt* + -| f"t - f + f"5$ +5 = 0 (9) 

This equation holds for arbitrary values of £. Equation (9) cannot be solved 
directly because of the unknown 5~derivatives. To obtain approximate values for 
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$£ and $>£, an approximate form of equation (9) is used in which only the terms 
involving £ and the derivative of the highest order are retained (in accordance 
with the theory of differential equations with large parameters) (ref. 4): 

5 = 0 ( 10 ) 

for which an analytical solution exists (ref. 4): 

*'( 5 , n ) - 1 - e "^* 11 ( 11 ) 

*(C,n) = n + — (12) 


The derivatives of equations (11) and (12) with respect to £ are 


and 



(13) 


(14) 


The derivatives from equations (13) and (14) are now inserted into equation (8) , 
which can then be solved by using the same numerical technique applied for the solu- 
tion of each equation of system (8) . The results are assumed to be valid not only 
for extremely high values, but also for moderately high values of 5 . Results of 
these calculations will also be discussed in direct comparison with the corresponding 
numerical results of the finite-difference procedure. 


NUMERICAL CALCULATION PROCEDURE 


The numerical procedure follows the ideas of Hall (ref. 10), who developed a 
method to solve the unsteady 2-D boundary layers for a flat plate impulsively set 
into motion. For the present case of harmonic oscillations, the time parameter, T, 
is first referred to the wavelength of the harmonic oscillation 


T = 


2tt/o)* 


(15) 


which changes equation (2) into 


3u 

3T 2lT 


+ u 


3u . - 3u 


3U u* 3fu 

3T 2it 3n 2 


(16) 


For discretization of equations (1) and (16) , a Crank-Nicolson type of scheme is 
used (fig. 2): information at points Q) , (2), and (3) is assumed to be given, and 

new values are calculated at position (?) . 
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This method is implicit in the direction normal to the wall (redirection > 
index n) and explicit in the two flow directions (T, index m and x, index SL ) . 

For central differentiation about point (5) in figure 2, the following expres- 
sions are derived with the abbreviations: 

. U n ~ 4 ^ U £+i,m,n + U £,m+i,n + U £,m,n^ 


l 

u = u„ , - u„ . - u„ 

n £+i,m,n £,m+i,n £,m,n 


m 

n ~ £,m+i,n £+i,m,n £,m,n 


(17) 


The velocities and their derivatives at midpoint (5) are 

_ 1 . s 

u g,+ (i/ 2 ) ,n+(i/ 2 ) ,n' 4 u £+i,nri-i,n U n 

(du \ _ _1_ , JL 

' 3x 4+(l/*),nrt-(l/2),n " 24x U ^ +1 ,m+1 ,n u n 

/9u\ _ 1 fl _1 s _ s \ 

' l9n '£+( 1 / 2 ),mf(i/ 2 ),n 2An H 1 » m+1 » n+1 " 4 U S-+i ,m+i ,n-i U n+i “ U n-i/ 

/Vu\ _ 1 /I • ^ 1 1 

\3n 2 A + (i/ 2 ),m+(i/2),n Ari 2 \ 4 U £+i ,m+i ,n+i 2 U £+i,m+i,n 4 U £+i,m+i,n-i 

, s 0 s s \ 

+ u , - 2u + u I 

n+i n n-i / 

( 5 ). 


^ u £+i,m+i,n + U n^ 


£+(i/z) ,m+(i/ 2 ) ,n 


(18) 


The introduction of the various terms of equation (18) into the momentum equation (16) 
yields 


a + b u +cu =d 

n £+i,m+i,n+i n £+i,m+i,n n £+i,m+i,n-i n 


(19) 


with the coefficients 


a = 


1 \ 


n ~ 8 An £+( 1 / 2 ) ,m+(i/ 2 ) ,n 4A n 2 


b = 


K , 1 

+ — U 


n 2AT 4Ax A+1 » m » n 2An 2 


c = 


n “ 8 An £+( 1 / 2 ) ,m+(i/ 2 ) ,n 4^2 j 


( 20 ) 
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\ ( 20 ) 
(cont) 


, - — m 1 2 _ 1 s £ _ 1 / s _ S .1 

n ~ 2^5 Ur 8Ax U ^+ 1 »®+ 1 »n 2 Ax n U n 2 Ap £+(1/2) ,m+(i/2) ,n ' n+i n-i' 1 

+ ^ (Vl ,*fi + + 2^ (I V. >m *l + “I <°*+> ,«fl + "*> 

, 1 / S „ s s . 

An 2 n+1 n n_1 

and with K = u)*/ 2 ir. 

For the continuity equation, the following expressions are derived: 

(— ) =i R— ) + (— ) ]) 

' 9x '£+(l/2) ,m+(i/2) ,n+(i/2) ^ L' 9 x A+(i/ 2) ,m+(i/2) ,n+i ' 9 x '£+(i/ 2) ,m+(i/2) ,nJ ( 

/ 9 V\ 1 r -I 

^ n 4 + (i/2),m+(i/2),n+(i/2) = An V ^(i/2),m+(i/2),n +1 " V l+d/2) .^(1/2) ,n 

An j. 

V £+(i/ 2) ,m+(i/2) ,n+i V ilH- (x / 2 ) »m+(i/2) »n 4 Ax U £+i ,m+i ,n+i 

. , £ , £, 

+ u„ , . +U.+UJ 

£+i,m+i,n n+i n 


( 21 ) 


To start the numerical calculation, simple, linearly extrapolated values from 
the positions (T) , (2) , and (5) are determined in position (4) . Equation (21) is 
first calculated for the normal velocity component through the boundary layer. With 
these initial estimates for v, equation (19) combined with the coefficients from 
equation (20) can be calculated by means of the solution of a linear system of equa- 
tions with a tridiagonal coefficient matrix (ref. 8). 

The newly calculated n values then enter equation (21), thus starting 

a second. cycle of the iterative procedure. This process is repeated until the dif- 
ferences in the velocity components, u, between two cycles are less than a given 
small number. Calculation then continues in the next grid point. 


The calculation first progresses in the x-direction (index £) until the maxi- 
mum x-value is reached and then continues one step, AT, in the time direction 
(index m) , etc. For the grid points (£, 1) at T « 0, the Blasius profiles, f’, are 
used as initial values. At the grid points (1, m) at x = Ax, the solution of equa- 
tion (8) is used as initial values along the time axis. 


Calculations were continued up to the third cycle, although the differences 
between the second and third cycle turned out to be negligibly small in most cases. 
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DISCUSSION OF THE RESULTS 


A series of calculations have been performed to test the validity and effective- 
ness of the numerical procedure. Comparisons with the corresponding low- and high- 
frequency solutions were made to investigate the accuracy of the results. 

Figures 3 and A show, first, wall shear-stress (x w ) and boundary-layer- 
displacement thickness (5 x ) over a period of time with 5 = x • gj* (x = the coordi- 
nate from the plate's leading edge) as a parameter. The numerical results for the 
first and second periods are compared with the low-frequency approximation. The 
amplitude of oscillation is e = 0.2. Except for some minor deviations at x « 0.1, 
agreement between both results is very good. Both t w - and 6 ^distributions show the 
well-known phase lead with increasing value of 5 (= xw*) , which tends to 45° for t w 
and gives a maximum value for up to approximately 20°. Figures 5 and 6 show 

corresponding results for x w and 6! plotted versus 5 for various dimensionless 
time levels. The differences between the first and second period of the numerical 
calculation decrease as time increases. Again, agreement between numerical results 
and the low-frequency solution is very good. 

Figures 7 and 8 show corresponding results for the higher-frequency case w* = 5. 
In figure 7 at approximately x = 0.6 (5=3) the wall shear stress has negative 
values over part of the cycle. The region of negative x, that is, the back-flow 
region, increases as x increases. The numerical calculation method proceeds 
through these reversed-flow regions without any numerical problems. For the largest 
x-values (x = 0.9), the results of the high-frequency approximation are discussed. 

Only minor differences in both cases are seen for x w and for 6^ In figure 8 (<S X ) 
results are plotted for the first three periods of oscillation. The number of cycles 
must be increased as the distance from the plate's leading edge increases to reach 
the steady, periodic solution. At 5 = 4.5 as many as four cycles must be calcu- 
lated until a good fit with the high-frequency results is obtained. 

Numerical calculation has been continued in the x-direction as far as x = 2.5, 
corresponding to a frequency parameter of 5 = 12.5. The calculation ran without any 
complications, although the back-flow region was extended to 0.4 < T < 0.8. However, 
with increasing back-flow velocities, the numerical calculation had an increased ten- 
dency to oscillate. These oscillations could easily be avoided by applying a simple, 
linear- interpolation procedure over three mesh points in the x-direction, which has 
already been applied successfully for calculation of 3-D steady boundary layers 
(ref9. 8, 9). \ 

There should be no problems for linear stability in cases in which the mesh 
sizes in the x- and T-directions are equal. The Courant-Friedrich-Levy (CFL) con- 
dition leads to the limit 


M > (22) 

with min(u£ +1 m n ) as the largest negative u-velocity component within the back- 
flow region, it’ is obvious that for Ax = AT: 

l min(u »+l,m,n ) l : 1 (23) 


The limit for equation (23) was by no means exceeded in the previously discussed cases. 
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Figures 9-12 show results for the high-frequency case (u* = 5) with increased 
oscillation amplitude e. It can be observed in the lower-amplitude case (e = 0.2, 
fig. 8) that the time dependency of deviates with an increasing value for x 

from a simple sinusoidal behavior: beyond the reversed-f low region the boundary 

layer steepens. This behavior of 6 X is much more pronounced with increasing oscil- 
lation amplitude and can be studied in detail in figure 10. If a Fourier analysis of 
S-lCT) is carried out, strong in-phase components of the second harmonic are obtained. 
The wall shear-stress curves (fig. 9), however, show only small deviations from a 
pure sinusoidal behavior. 

Figure 11 shows details of the u-velocity components in the vicinity of the 
wall where considerable back flow occurs over part of the cycle. Correspondence 
between the numerical data and the high-frequency results are, aghin, extremely good. 
This agreement ensures that the numerical method is of sufficiently high accuracy and 
reliability. 

To investigate the influences of increasing oscillation amplitude, e was 
increased to values as large as 0.5, 0.7, and 0.9. Figures 12 and 13 show the 
corresponding high-frequency results for t w and 6^ It is remarkable to see the 
increasing peak of beyond the back-flow region in figure 13. At e = 0.9 the 

peak value reaches 14.6, and at e = 1.0 the calculation discloses a singularity at 
point T = 0.75. This singularity is caused by the definition of the displacement 
thickness which refers to the velocity at the outer edge of the boundary layer. This 
velocity tends to zero at a particular point of the cycle for the case of e = 1. 

The wall shear-stress behavior (fig. 12) is regular as in the smaller e-cases, with- 
out irregularities even at e = 1. 

The numerical-calculation procedure, however, was limited as soon as the 
reversed-flow region reached back to the leading edge of the plate during part of the 
cycle. This was the case for to* = 5 at e = 0.5 and 0.65 > x > 0.75. The solu- 
tion then started to oscillate with increasing amplitudes. There was no possibility 
of avoiding these oscillations by means of an interpolation procedure. It is assumed, 
however, that the breakdown of the numerical method is caused by difficulties at the 
leading edge of the plate, and not by the increasing, reversed-flow regions. 


CONCLUSION 


A numerical method has been developed to calculate the unsteady laminar boundary 
layer over an oscillating flat plate. This investigation was made to gain insight 
into flow behavior, and to establish the applicability of a numerical-calculation 
procedure which, based on the boundary-layer concept, could be applied as soon as 
flow reversal occurred over parts of the cycle. Corresponding low- and high-frequency 
solutions have been developed for comparison with the full numerical results. For 
both of these limiting cases, agreement with the numerical results were good. No 
numerical problems occurred within regions of reversed flow. The numerical calcula- 
tion procedure, however, came to a limit as soon as the back-flow region reached the 
leading edge of the plate. It is assumed that in this case the boundary layer con- 
cept is violated because of the more complicated flow behavior at the plate's leading 
edge. The most remarkable behavior is shown by the boundary-layer displacement thick- 
ness over a cycle of oscillation. With increasing amplitude, increases rapidly 
beyond the back-flow region. This behavior can be observed with the relatively simple 
high-frequency approximation discussed previously. 
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The numerical procedure developed for oscillating flat plates should permit the 
development of a similar method for using oscillating profiles to investigate unsteady 
separation at high incidences and at moderate to large oscillation amplitudes. 
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Figure 2.- Finite difference grid. 
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Figure 3.- Development of wall shear stress with respect to time at low frequency 

(to* = 1, e = 0.2). 
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Figure 4.- Development of boundary layer displacement thickness with respect to 

time at low frequency (to* = 1 , e = 0.2). 
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Figure 7.- Wall shear stress with respect to time at high frequency (w* = 5, e = 0.2). 
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Figure 9.- Wall shear stress with respect to time for to* = 5, e ~ 0.3. 
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Figure 10.- Boundary layer displacement thickness with respect to time for 

u* = 5, e = 0.3. 
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Figure 12.- Influence of oscillation amplitude on wall shear stress (w* = 5, x = 1.1) 
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Figure 13.- Influence. of oscillation amplitude on boundary layer displacement 

thickness (oj* = 5, x = 1.1). 
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Calculations, however, came to a limit where the back-flow region reached 
the plate's leading edge in the case of high oscillation amplitudes. It is 
assumed that this limit is caused by the special behavior of the flow at the 
plate's leading edge where the boundary layer equations are not valid. 
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